#tables
#start with wksp_2010_09_03.rdata

library(biomaRt)
chick.desc <- getBM(attributes=c("ensembl_transcript_id","description"),filters="ensembl_transcript_id",values=grouse.chick.map$ensembl_transcript_id,mart=chicken.mart) 

#create table for Fst hits
#ssh_1675-1 missing from BBlikTest.pooled.df2, blast and find hit later
#or just ignore
#contig10 and contig62 are the same
write.table(merge(BBlikTest.pooled.df2[BBlikTest.pooled.df2$bonferroni<0.05,],chick.desc),sep="\t",quote=FALSE,row.names=FALSE,file="../figures/fst_table.txt")

#create table for top omega
for(icol in c("Dn","Ds","transitions", "transversions","dN","dN.se","dS","dS.se","kappa","omega","Ls","Ln","t" )){
	grouse.chick.MK[grep('Error',grouse.chick.MK$Ls),icol] <- "-1"
	grouse.chick.MK[,icol] <- as.numeric(grouse.chick.MK[,icol])
	}

grouse.chick.MK <- grouse.chick.MK[!is.na(grouse.chick.MK$grouse.id),]

#produces bizarre NAs
tmp <- grouse.chick.MK[grouse.chick.MK$omega>quantile(as.numeric(grouse.chick.MK$omega[as.numeric(grouse.chick.MK$omega)>0 & as.numeric(grouse.chick.MK$omega)<10 &!is.na(grouse.chick.MK$omega)]),0.95,na.rm=TRUE),]
#selects top 5%
tmp <- tmp[-grep("NA",row.names(tmp)),]
tmp <- tmp[tmp$omega < 10,]
write.table(tmp,file="../figures/omega_table.txt",sep="\t",quote=FALSE,row.names=FALSE)

#Tajima D tables, may not be the best way to select genes
#1% cut-off top > 3.58, 1% bottom < -1.57
write.table(grouse.TjD.df[grouse.TjD.df$pop1.max>quantile(grouse.TjD.df$pop1.max,.99)|grouse.TjD.df$pop2.max>quantile(grouse.TjD.df$pop2.max,.99)|grouse.TjD.df$pop3.max>quantile(grouse.TjD.df$pop3.max,.99),],file="../figures/TjD_top_table.txt",sep="\t",quote=FALSE,row.names=FALSE)

write.table(grouse.TjD.df[grouse.TjD.df$pop1.min<quantile(grouse.TjD.df$pop1.min,.01)|grouse.TjD.df$pop2.min<quantile(grouse.TjD.df$pop2.min,.01)|grouse.TjD.df$pop3.min<quantile(grouse.TjD.df$pop3.min,.01),],file="../figures/TjD_bottom_table.txt",sep="\t",quote=FALSE,row.names=FALSE)


#below added to Bayes_Fst, Bayes Fst also altered a little to test a greater range of Fst
BBlikTest.pooled.df$fst <- 0

BBlikTest.pooled.df$fst <- apply(BBlikTest.pooled[,-1],1,function(X){
	fst.seq <- 1/(seq(0.01,20,length=100)+1)
	fst.seq[which.max(X)]
	})
